Specific heat and high-temperature series of lattice models: interpolation scheme and 
examples on quantum spin systems in one and two dimensions 
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We have developed a new method for evaluating the specific heat of lattice spin systems. It is 
based on the knowledge of high-temperature series expansions, the total entropy of the system and 
the low-temperature expected behavior of the specific heat as well as the ground-state energy. By 
the choice of an appropriate variable (entropy as a function of energy), a stable interpolation scheme 
between low and high temperature is performed. Contrary to previous methods, the constraint that 
the total entropy is log(25+ 1) for a spin S on each site is automatically satisfied. We present some 
applications to quantum spin models on one- and two- dimensional lattices. Remarkably, in most 
cases, a good accuracy is obtained down to zero temperature. 

PACS numbers: 75.10.jm 



I. INTRODUCTION 

The accurate knowledge of the thermodynamic quan- 
tities of quantum magnets is an important issue from the 
experimental point of view. It allows one to determine 
precisely exchange energies J from experimental data or 
to identify possible deviations from a given model. In 
this paper we propose a simple method to compute the 
heat capacity from a high-temperature (HT) expansion. 
Numerous techniques have been developed to study ther*. 
modynamical singularities from power series expansions.!!! 
However, in this work, the point of view is slightly dif- 
ferent. We mainly focus on the computation of the heat 
capacity for systems without a phase transition at T > 
(here in one and two dimensions) and try to evaluate 
the specific heat c v (T) accurately in the largest possible 
temperature range. We devised a two-point Pade-like in- 
terpolation of a particular function, namely the entropy 
as a function of energy, which satisfies the energy and en- 
tropy sum rules obeyed by c v (T). These two constraints 
are non-local in temperature and improve significantly 
the convergence of standard Pade approximations. 

Before describing the method itself, let us briefly re- 
view some commonly used methods to investigate the 
thermodynamic properties of spin systems. 

High-temperature series expansion is the usual 
approach to evaluate the strength of microscopic inter- 
actions from experiments. Thermodynamic quantities 
such as magnetic susceptibility \ or specific heat c v are 
expanded in powers of = X/T. Presently computers 
do not allow to compute more than about 20 terms for 
quantum systems.El When analyzed through (differential) 
Pade approximants, HT series give reliable results for T 
greater than the typical magnetic exchange energy of the 
model. When the low-temperature form of a given quan- 
tity is known, Pade approximants can be improved to re- 
produce the correct asymptotic behavior. However, with 



this approach it is difficult to investigate temperatures 
much below J. In some cases, a detailed understand- 
ing of the low-energy regime allows one to construct bi- 
ased approximants and to compute thermodynamics in a 
larger temperature range (see for instance Ref. || for the 
magnetic susceptibility of square- and triangular-lattice 
ferromagnets) . Another series approach to thermody- 
namics is the finite-cluster-expansion. This expansion in 
powers of the coupling constant may, in some cases, al- 
lows one to go tpi lower temperature. It has been applied 
in onecl and twcQ'B dimensions. 

Finite-size calculation. The previous approaches 
are exact at high temperature but do not satisfy the sum 
rule that the total entropy is ln(25 + 1). An alterna- 
tive is to compute the heat capacity of a finite system 
from the exact spectrum of a small cluster of spins, and 
then compute the partition functionJiy summing over all 
eigenstates (example: spin-1 chainLrtl). A direct finita- 
temperature lanczos algorithm can also be performed!!, 
as well as transfer matrix techniques. These techniques 
automatically give the correct entropy, but they are lim- 
ited to small systems (N < 36) and finite-size errors may 
be difficult to control, especially in two dimensions. 

Quantitative extrapolations to N — oc of finite-size 
c v ,n(T) data have been done in spin chains but are not 
really efEcient.- at .l ow temperatures. It has been applied 
to spin chainsQE3EJIi3 and the triangular antiferromag- 
net.t^l Power-law behaviors at low temperature cannot be 
observed due to the important discretization of the low 
energy spectrum in a small system. 

Quantum Monte- Carlo simulations can reach a 
rather large number at spins for not frustrated systems 
(as the square latticeEfl). These calculations reproduce 
the correct HT behavior (up to small statistical errors) 
but cannot reach the very low temperature regime. This 
technique is one of the most efficient in studying thermo- 
dynamics properties but it requires important numerical 



2 



effort and is limited to small systems when the model 
contains frustration. 

Sum rules. We propose a simple method to compute 
the heat capacity which involves a new kind of HT series 
analysis. In contrast to previous methods, we provide a 
procedure for c v (T) which satisfies a first sum rule: 



c v (T) 
T 



dT = ]n(2S+ 1). 



(1) 



Moreover, we incorporate the knowledge of the ground 
state energy eo and the energy e m at infinite temper- 
ature. These energies are either known exactly or can 
be computed numerically. We use this information in a 
second sum rule: 



ln(2) 




FIG. 1: Typical shape of the entropy 
e for a spin-j system. 



s as 



a function of energy 



c v {T)dT = e(T = oo) - e(T = 0) 



(2) 



Finally, the low-temperature leading contribution of 
c v (T) is taken into account. 

We claim that the implementation of these integral 
constraints increases irua significant way the range of va- 
lidity of the HT series.t3 We will show that this method is 
indeed successful even at low or zero temperature for the 
models we have investigated. By "successful", we mean 
that the specific heat is obtained down to zero tempera- 
ture with a relative accuracy typically between 1% and 
0.1% with only ten terms in the high-temperature se- 
ries. This accuracy should be considered as surprisingly 
high for low-temperature quantities obtained with a high- 
temperature series expansion. 



II. INTERPOLATION PROCEDURE 

A. Elementary thermodynamics: From c„(T) to s(e) 

In the canonical ensemble, the energy per site e and the 
entropy per site s are increasing functions of temperature 
T = 4. Since e is a monotonic function of T, we can 
express s as a function of e rather than of temperature. 
The specific heat as a function of temperature is easily 
expressed with s(e) in the following way. The entropy 
s(T) obeys: 



ds c v 
dT ~ f 



(3) 



Where c v is the specific heat per site. It follows that : 



. ds c v dT 
de T de 



(4) 



where the prime denotes differentiation with respect to 
e. Since ^ = we have (the Boltzmann constant ks 
is set to unity in the whole paper) : 



s'(e) 



1 

T 



(5) 



We use this relation to eliminate the variable temperature 
in c v : 



de 


'dT' 


-l 


\-(-)\ 


dT = 


de 




de \s'J 



We eventually find 



(G) 



(7) 



Eqs. (J5|,[7|) are the basic relations that we shall use here. 
These relations can also be obtained by considering the 
density of states of a large but finite systemJl3. For sim- 
plicity, we^discuss the case of a quantum spin-i model on 
a latticeO The simplest cases are the Heisenberg models: 



H 



(8) 



<hj> 



where the sum runs over the first nearest neighbors. s(e) 
is defined inside an interval going from the ground-state 
energy eo up to the high-temperature energy e m . e m 
is the free-spin average of the energy, which is obtained 
straightforwardly since, at infinite temperature one has 

Si ■ Sj^ = 0. For the Hamiltonian (||), one simply gets 
e m = 0. If the ground state |0) is ferromagnetic ( J < 
in Eq. (|)), one has (0| S, ■ §j |0) = ± and e = Jf for 
a Bravais lattice of coordination number z. However, 
for antiferromagnetic models, eo is not known exactly 
but Monte-Carlo simulations, exact diagonalizations or 
analytical calculation (spin-wave, mean-field Schwinger- 
bosons, etc) can be used. In summary the function s(e), 
defined in the interval [eo,0] (eo < 0), is an increasing 
function of e, starting at in eo with an infinite slope, 
and finishing at ln(2) in e = (see Fig. |l|). 

To find an approximation of s(e), we combine three 
types of exact informations on s(e): 

a) The entropy per spin is ln(2) at infinite tempera- 
ture 



s(e m = 0) = ln(2) 



(9) 
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and vanishes at zero temperatureS 
a(eo) =0 



(10) 



b) The behavior of c v (T — > oo) is known from HT 
series. From the expansion of c v (T) in powers of 
1/T up to order l/T" , we get the expansion of 
s(e) in the vicinity of e — ► e m = in powers of e to 
order e": 



s(e) f 



ln(2) + ^a. ( e 4 



(11) 



i=2 



This expansion can be computed by solving 
Eqs. (^,0) order by order for s(e) (details in ap- 
pendix A). 

c) When the low-energy physics of the model is under- 
stood, the low-temperature limit of the specific heat 
can often be predicted. In the case where the sys- 
tem exhibits some ferro- or antiferromagnetic long- 
range order at zero temperature, the low-lying ex- 
citations are spin waves. These gapless modes give 
a low-temperature heat capacity which is a power 
law. When the space dimension is D, a ferromag- 
net has c v (T) ~ T D / 2 and an antiferromagnet has 
c v (T) ~ T D . In both a cases, we can write: 



c v {T) 



(12) 



where p and q are integers. This low temperature 
behavior translates into a behavior of s(e) about 
the ground-state energy per site eo. 



s(e) e _> eo - (e - e )p+" 



(13) 



On the other hand, if elementary excitations are 
gapped the system has a thermally activated spe- 
cific heat 



c v {T) 



exp(-A/T) 



(14) 



and Eq. ( p^ ) is replaced by a logarithmic behavior 
about the ground-state energy (see appendix C). 



B. Interpolation by Pade approximants 

We now define the interpolation procedure between low 
and high energy for s(e). We look for an approximation 
of s(e) which satisfies Eqs. (|), ©, (|ll|) and (|l|). A 
two-point Pade interpolation is not directly possible since 
s(e) is singular in e = eo- We have to transform s to 
obtain regular behavior in the whole interval [eo, 0]. This 
is possible when p and q are two integers (but other low- 
temperature form will be used for thermally activated 
c.„ in gapped systems, see section [V). We define a new 
function G(e) as: 



G(e) = s(e) 



p+q 



(15) 



This function now behaves as (e — eo) p when e — > eo 
and as ln(2) p+9 + 0(e 2 ) when e — > 0. It is now possible 
to look for approximations G app of G which are of Pade 
form (details in appendix B). This interpolation scheme 
fails if any of the functions G(e), s(e), s'(e) or — s"(e) 
becomes negative in the interval [eo , 0] . This provides 
a natural criterion to select the degrees [u, d] for which 
Pade approximants give physical solutions. The curve 
for the specific heat c v (T) is eventually obtained in a 
parametric form {T(e), c 1 ,(e)} ee [ eo Q i thanks to Eqs. (Q) 
and (HI). 



III. LOW-TEMPERATURE POWER-LAW 
BEHAVIOR OF c v (T -> 0) 

We start by illustrating the method for systems where 
the specific heat is proportional to T a at low temper- 
ature. Some cases where the specific heat is thermally 
activated will be described in the next section. Data files 
with the numerical results presented here (as well as col- 
ored versions of the figures) are available upon request. 



A. One-dimensional S = i XY model 



As a first application, we consider the spin-^ XY chain: 



H = 2 



[pi ' ^ 



i+1 



oV qV 



(16) 



This spin problem can be solved exactly (the Jordan- 
Wigner transformation maps the XY model to free spin- 
less fermions) and provides a check for a gapless spec- 
trum. The energy per site at T = (resp. at T — > oo) 



IS Gq. 

byL s 



-- (resp. e m = 0). The heat capacity is given 



Cv(J3) = ^ 

7T 



tt/2 



cos 2 (fc) 



o cosh (/3cos(fc)) 



-dk 



(17) 



From this formula, the high- and low-temperature ex- 
pansions can be computed. We find: c v ([3 — ► 0) = 
i/3 2 - |/3 4 + • • • and c v (T -> 0) = f T + 0(T 3 ). The lin- 
ear law at low temperature gives p = q = 1 in Eq. (|T^). 
Fig. U shows the comparison between the exact c v (T) and 
the ones obtained from the Pade approximants. At each 
order n, a few approximants [it, n — u] lead to the same 
variations. In Fig. ^, only diagonal Pade approximants 
[§, j] are shown. A convergence of the specific heat is 
obtained with a relative accuracy of the order of one per- 
cent down to zero temperature with only six terms in the 
HT expansion. We want to emphasize that the slope 
^f(T = 0) is not imposed in our procedure. There- 
fore, the method provides quantitative information on 
this low-energy parameter which characterizes the low- 
energy excitations. The value of the prefactor is exactly 
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0.1 1 
T 



FIG. 2: One-dimensional S = § XY model. Top: Ex- 
act specific heat (full line), approximations at orders n = 2 
and n = 4 (dash-dot lines), n = 6 (dotted line). Bottom: 
Differences with the exact result (notice the magnified scale), 
n = 6 (Dotted line), n = 8 (dash line), n = 10 (long dash), 
n = 22 (dash-dot), n = 30 (full line). All Pade approximants 
are diagonal: [n/2, n/2]. 

known for this XY model and we find that it oscillates 
around the exact value (i.e. 7r/6). It is a few percent 
below the exact value at the order n = 6. 



B. One-dimensional antiferromagnetic Heisenberg 
spin-| model 



H = 2Y / S i -S l+1 (18) 

i 

This model is solvable with the Beth^ ansatzS The 
ground-state energy is exactly knowncil to be eo = 
— 21n(2) The low-temperature limit of the heat- 

capacity isE3 c v (T) = | -tjp(T 2 ). The HT expansion 
computed by Baker et a/.cp-,was recently extended to 
order /3 24 by Biihler et al.E-i The heat capacity curve 
has been computed for anisotropic versions of the spin- ^ 
chained, but we are not aware of any exact computation of 
c v for the isotropic Heisenberg model based on the exact 
Bethe ansatz equations (due to the existence of an infi- 
nite number of coupled equations at the isotropic point). 
However, an .approximate solution was proposed recently 
by Klumper.EEl The accuracy of specific heat obtained by 
his method is claimed to be extremely high (~ 10~ 7 ). 

Our results for this model are displayed Fig. ||. A good 
convergence is obtained down to zero temperature. If we 




0.01 0.1 1 10 

T 



FIG. 3: Specific heat of the antiferromagnetic soin-^ 
Heisenberg chain. Top: Result from Klumper et a£3 (full 
line), approximations at orders n = 3 and n = 5 (dash-dotted 
lines) and n = 7 (dots). Bottom : Difference between RefE3 
and present approach, n = 7 (dot line), n = 9 (dash line), 
n — 11 (long dash), odd n from n = 13 to n = 23 and n — 24 
(full lines). 



truncate the HT expansion at order seven, an absolute 
precision of one percent can already be obtained. At the 
highest order, the specific heat is in excellent agreement 
with the result of Ref. ^6|. The height of the peak, for 
instance, agrees with an accuracy of 1CT 4 . The low tem- 
perature limit is also in good agreement with the exact 
result: we find c v (T) ~ 0.329T at order 24 (instead of 
Cv (T) = §T). The difference with result of Ref. |^ at in- 
termediate temperatures (see bottom of Fig. ||) is due to 
the fact our the low-temperature coefficient is off by 1% 
from the exact value ^. This result is obtained without 
any extrapolation when n — * oo. Taking into account 
finite-order corrections could improve further the accu- 
racy. 



C. Triangular lattice Heisenberg models 

First, we look at the convergence of this method on 
the ferromagnetic case where the ground state energy is 
known exactly. Fig. ^| shows the comparison of the vari- 
ous Pade interpDlants [u, d], with n = u + d = 13 (highest 
available ordercH). Two of them have a singularity ([10, 3] 
and [4, 9]). If we exclude the polynomial form (d = 0) and 
the cases d > 10, all interpolations lead to the same vari- 
ations of c v (T) in the whole range of temperature. Fig. ^ 
also illustrates the very small dispersion of the differ- 
ent Pade approximants obtained with this method. It is 
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FIG. 4: Ferromagnetic spin-1/2 Heisenberg model 
on the triangular lattice. Comparison the specific heats 
obtained from different Pade approximants at order n = 13. 
The Pade approximant of reference has degrees [7, 6]. Approx- 
imants [6, 7], [8, 5], [12, 1] and [3, 10] differ from the reference 
by less than 10" 3 . 

also remarkable that the coefficient of the dominant term 
at low temperature depends weakly on the interpolation 
function used. As before, only the power of T is imposed, 
but the prefactor is not. We get c v {T) ~ .142(2)T at low 
temperature. Notice that this quantity is not given cor- 
rectly by a linear spin-wave approximation (non interact- 
ing magnons). We are not aware of any previous result 
for this quantity. A high-density mono-layer of solid 3 He 
solid is, to our knowledge, the only experimental realiza- 
tion of a triangular-lattice spin-i ferromagneL The heat 
capacity has been measured by Ishida et. a£Q and their 



results are in very good agreement with our calculations. 

In Fig. |5|.a and Fig. [jj.b, we see now how the result 
converges when more terms are added in the HT expan- 
sion. For each HT expansion order n (as indicated in 
the figure), all possible fractions [u, d] (u + d = n) are 
tried. Except for the lowest order, at each order several 
approximants fall on the "same" curve. A reasonable 
convergence is obtained even for n = 5. 

Fig. g.b shows the approximated c v (T) with increasing 
HT-order n for the antiferromagnet Hamiltonian. jEfois, 
model is frustrated, but it is now well established 2 ^'!^ 
that the ground-state is a three-sublattice Neel state. 
The low temperature specific heat is thus proportional to 
T 2 . The ground-state,energy was estimated from exact 
diagonalization data.Ej Again here, we see a convergence 
for n > 5. This should be compared with direct Pade ap- 
proximants to the specific heat, which was done in Ref. |27] 
(dotted curve in Fig. |b). The comparison with our 
results shows indeed that the direct approach does not 
even allow to get correctly the maximum of c v (T) . Our 
method converges with a relative accuracy of about 1% 
down to zero temperature and we get c v (T) ~ 5.3(2)T 2 . 
It would be interesting to compare this result with spin- 
wave calculations taking magnon-magon interactions into 
account. 



D. Square lattice Heisenberg model 

We evaluated the specific heat for the S = h Heisen- 
berg model on the square lattice. Fig. || shows the con- 
vergence of the specific heat with respect to the number 
of terms in the HT-expansion. The convergence with 
the HT-expansion order is faster in the antiferromag- 
netic case than in the ferromagnetic one. At low tem- 
perature we obtain c v (T) = 0.25(0.01)T (resp. c v (T) = 
0.25(0.01)T 2 ) for the ferromagnetic (resp. antiferromag- 
netic) case. The result in the antiferromagnetic case can 
be coppared with the modified spin- wave theory of Taka- 
hashM which is c v (T) = 0.214T 2 in our units. 



IV. SYSTEMS WITH THERMALLY 
ACTIVATED HEAT CAPACITY 

A. One-dimensional Ising model 



H = 2j2Sf-S? +1 (19) 

i 

This model is the simplest one with a gapped spec- 
trum. It allows to check the convergence of our method 
to the exact result in the case of a low temperature ac- 
tivated heat capacity. The energy per site at T = is 
eo = —1 and e m — at T = oo. The specific heat per 



spin is given by c v {(3 = 1/T) = (/?/ cosh(/3)) 2 . The HT- 
expansion starts as c v (f3 = 1/T) ~ /3 2 — /3 4 + . . ., while 
the spectrum has a gap A = 2. The heat capacity is 
thermally activated at low temperature: 

c v (/3 -> oo) ~ (3 2 exp(- A/3) (20) 

For such behavior, the entropy behaves as 

s(e) ~ -(e-eo)log(e-e )/A (21) 

when e — > eo- We need a transformation which converts 
this logarithmic singularity into a regular behavior. We 
choose the transformation 

G(e) = (e-e )±(^-) (22) 
de \ e — eo / 
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FIG. 5: Spin-1/2 Heisenberg model on the triangular lattice, (a) - Top: Variations of c v (T) with respect to the number 
of terms in the HT expansion in the ferromagnetic case. The even n from 2 to 8 and n = 13 are shown. Inset: data for 
n — 13 in a log-log plot. Bottom: Difference with the highest order (n = 13). Even n from 6 to 12 are shown. Order 12 and 13 
differ by less than 3.10 -3 . (b): Same as (a) for thejantiferromagnetic case. Order 12 and 13 differ by less than 2.10" 3 . The 
dotted line (with the spike) is from Elstner el. al£3 




FIG. 6: spin-^ Heisenberg model on the square lattice. Same as Fig. g for the (a) square lattice ferromagnet and (b) 
square lattice antiferromagnet. Even n from 2 to 14 and n — 13 are shown. Orders n — 13 and 14 differ by less than 4.10 -3 
in the ferromagnetic case and by less than 10 -3 in the antiferromagnetic case. 
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Model 


F/AF 


eo — e m 


n 


c(T->0) 




max 


Heisenberg 
chain S = i 


AF 


ff n 


— 1 1 

24H 

rn 


0.329T 


0.9618(2) 


0.3497 


Heisenberg 
chain S = 1 


AF 


-1.401H 




exp (-A/T) 
A = .40* 


0.861(1) 


0.543(1) 


Heisenberg 
triangular lat. 


F 


— iTT- 


130 


0.142(2)T 


1.375(5) 


0.403(3) 


AF 


-1.11S 


5.3(2)T 2 


0.84(1) 


0.2231(5) 


Heisenberg 
square lat. 


F 


-1 


140 


0.25(0.01)T 


0.785(4) 


0.383(3) 


AF 


SB 1 ' 34 


0.25(0.01)T 2 


1.163(2) 


0.467(2) 



TABLE I: Parameters of the different models. Column two indicates whether the couplings are ferromagnetic (F) or antifer- 
romagnetic (AF). eo is the ground-state energy per spin, e m is the mean energy at infinite temperature (here e m = 0). n is 
the highest known order of the HT expansion of the specific heat. The three last columns are parameters extracted from our 
analysis: low-temperature limit and position of the peak. These result are obtained from the highest order available. The error 
bars reflects the dispersion of the different Pade approximants at the highest order. Figures with a star (*) are those for which 
the exact value is known (see text). 



B. One-dimensional 5=1 Heisenberg model 




FIG. 7: One-dimensional Ising model. Top: comparison 
between the exact specific heat (full line) and those obtained 
from second (dot-dashed line) and third (dashed line) orders. 
Bottom: difference between exact result and the approxima- 
tions at orders n — 2, 3, 4 and 5. 



and we approximate G by a Pade approximant (see ap- 
pendix C for details). The convergence is shown in fig- 
ure [7| . The gap also converges rapidly to the exact value 
2. The best Pade approximants are those of the form 
[d,d] with d~ f. 



H 



5 Si ' Si 



+1 



(23) 



It is known that this system exhibits a spin gap (Hal- 
dane phase) and the heat-capacity is thermally activatad- 
More precisely, a non-linear sigma model approaches 
gives: 



c v {T) 



'2tz 



exp (-A/T) 



(24) 



More generally, if elementary excitations are weakly 
interacting massive bosons, one expects c v (T) ~ 
exp(— A/T)T-°/ 2 ~ 2 in space dimension D. 

The value of the ground-state energy and spin gap 
is known very accurately (from the density-matrix 
renoenaalization-group calculations done by White and 
HusdB) to be e = -1.401484038971(4) and A = 
0.41050(2). We use the HT series data obtained up to 
order 20 by Elstner, Jolicceur and Golinelli.a We ap- 
ply the same transformation as with the ID Ising model. 
Notice that Eq. @ implies c v (T) ~ exp (-A/T)/T 2 in- 
stead of- exp (-A/T)/Tt as it should be from Eq. (g). 
The consequence of this little inconsistency is that, at 
fixed order in the HT series, we overestimate the value of 
the gap A. Note that the correct low-T behavior would 
only give a logarithmic correction to s(e) (Eq. |T]). This 
can be rectified, in principle, by a more sophisticated 
transformation of s(e). Nevertheless, we obtain a gap 
which varies roughly as 1 /n 2 (n is the order of the HT- 
expansion) and the extrapolation shown Fig. || gives 0.40. 
So, very surprisingly, this HT approach is able to provide 
some quantitative and non-trivial information about the 
low-energy physics. 
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FIG. 8: Specific heat of the antiferromagnetic Heisen- 
berg spin-1 chain. Top: HT orders are n = 2, 4, 8 and 
20 (full line). Inset: Log- log plot of n — 20 and theoretical 
prediction of Eq. ( p4j ) (dashed). Bottom: Difference between 
n = 20 and even n from 8 to 18. n = 16, 18 differ from n — 20 
by less than 2.10" 3 . 



V. CONCLUSIONS 



FIG. 9: Convergence of c v (T) at low temperature for the an- 
tiferromagnetic Heisenberg spin-1 chain. Even orders 
from 8 to 20 are showg. Full bold line is n = 20, the dashed 
bold line is from RefB (Eq. Notice that n = 16 and 18 
are almost identical at the scale of the figure. Inset: Spin gap 
for orders n = 8, 10, 12, 14, 16, 18, 20. The horizontal dashed 
line is the exact value (0.41) and the full line is a least-square 
fit to the data. It extrapolates to A = .40. Upward triangles 
are for approximants with degrees [ J + 1, § — 1] , square for 
[f , f ] and downward triangles for degrees [2 + 



We have presented a new and simple method to ana- 
lyze HT series expansion for the specific heat of spin sys- 
tems. It requires only : the ground-state energy and the 
qualitative behavior of c v (T) at low temperature. These 
two pieces of information allow us to constrain the spe- 
cific heat with two sum rules which improve drastically 
the convergence of the HT expansion. This technique 
is particularly appropriate to analyze specific heat mea- 
surements in a large temperature range. Since in many 
cases it is able to predict the position and the height 
of the maximum of c v (T) , this method can simplify and 
improve the determination of exchange parameter(s) and 
of the number of spins in the sample. Remarkably, the 
method converges down to zero temperature with only 
ten terms of the HT series in most of the cases we have 
investigated. This method can even provide some quanti- 
tative information on the low-energy physics (e.g. of the 
value of the spin gap of the Haldane chain, or the zecp- 
temperature entropy of the Kagome antiferromagnetE3). 

It would be interesting to apply this technique to more 
sophisticated models such as the t— J or Hubbard models. 
It also seems worthwhile to investigate other interpola- 
tion methods than the simple Pade approximants pre- 
sented here. The application of this new HT analysis to 
systems with finite-temperature phase transitions is an- 
other promising direction we are currently investigating. 
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APPENDIX A: SERIES FOR s(e -» 0) 

We describe how the series for s(e — > 0) is obtained 
from the series for c v (T — * oo). Assume the expansion of 
the specific heat is known up to order n: 



c v {T)i 



E 

»=2 



(25) 



Eqs. 



and (I?]) implies 
s"{e)c v (T = l/«'(e)) 



s"(e) 



[s'(e)f 



(26) 



Expanding Eq. (|26|) in powers of e gives the series for 
s(e — > 0). For instance, if we have 

1 , 



c v (T)q 



«4 



2^2 2^3 ^4 



0{ 



T 5 ' 



(27) 



9 



we obtain 



S (e)_ =ln(2)--Le 2 --5e 3 - 
2a 2 



24a| 



3a| e 4 



In fact, solving Eq. (lit 
a software like Maple: 
4; Order 



(28) 

) can be done very simply with 



n := 4; Urder:=n+1: 
cv := add ( a[i]/T**i, i=2..n ); 
eq := -D(s)(e)**2 / D(D(s))(e) 

= subs(T=l/D(s) (e) ,cv) ; 
dsolve({eq, s(0)=ln(2), D(s)(0) = 0} 
s(e) , ' type=series ' ) ; 



Pade [M] (e) 



1 + - 

p 



— p — 14- — 



(31) 



Pade [0 . 2] (e)- 1 = 1 - p— + \ [p(p - 1) + x] ( - ) ( 82 ) 
e 2 Ve / 



where 



APPENDIX B: PADE APPROXIMANT FOR G(e) 

When the system has a specific heat with a power law 
T? at low temperature the function G(e) = s(e) p+q 
regular at the ground-state energy e " 
high temperature e 
Gby: 



is 



eo as well as at 
0. We can therefore approximate 



G app (e) = ln(2) 



p+q 



1-- Pade [M] (e) (29) 



where Pade^] (e) = N(e)/D(e), N (resp. D) is a poly- 
nomial of degree u (resp. d) and -/V(O) = -D(O) = 1. By 
construction, Eq. (p9|) guarantees that G app has the cor- 
rect behavior at low temperature G app ~ (e — eo) p . One 
has to evaluate the Pade approximant so that the HT ex- 
pansion of G(e) (and thus of s(e) and c„(T)) is exact up 
to order n = it + d. As usual, this is done by expanding 
both sides of Eq. ( |29| ) in powers of e and solving a linear 
system to determine the unknown coefficients of the two 
polynomials N and D. As an example, we give here the 
general expression for the three possible approximants at 
order n = 2: 

Pade [2 , 0] (e) = 1 + p- 4- ~ \p(p + 1) - x] (-) (30) 
e 2 \e J 



x = 



(33) 



a 2 ln(2) 

APPENDIX C: G(e) FOR GAPPED SYSTEMS 

With gapped systems, the specific heat has the form: 
c v (T) = exp(A/T)T a ' 2 . We obtain the low energy limit: 



s(e) 



eo 



, ,, \ * \ -, ,—ln(e — eo) , 
ln((e - e )A) + aln( y — ^) 



_ (34) 
For a = 0, we recover Eq. ( pi| ) and the transformation 
of Eq. ( |2^ ) holds. When a ^ 0, in order to remove the 
logarithmic singularities, one has to differentiate several 
times. But such transformations provide singular behav- 
ior for s(e), mainly because of the change of sign in the 
arguments of the log between e = eo and e = e m = 0. 
More work has to be done in these cases. 

Unlike the case of a power-low behavior at low temper- 
ature, cases with a gap require an integration in order to 
go back from G(e) defined Eq. ( p2| ) to the entropy and 
the specific heat. Since we look for approximations of 
G(e) in a Pade form, this integration can be performed 
analytically. The value of the gap A can be obtained di- 
rectly from G, without integration, since G(e = eo) = t. 
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